/***********************************************************************
Compile       : gcc -O3 -o checktraj1 checktraj1.c dynamics.c ../design/design.c ../design/utils.c -I../design
Usage         : ./checktraj1 mode mode1 mode2 irflag N L inputfile 
                (./checktraj1 0 0 0 1 11 12 Budding.txt)
mode          : 0 normal ; 1 cycled ; 2 triggled cycled
mode1         : 0 disallow self-activation ; 1 allow self_activation
mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II
irflag        : 1 irreducible solution only
N             : number of the nodes
L             : number of the steps        
Output        : activity min_link designability1 designability2
activity      : number of 1 in the process
min_link      : minimum links of the qualified network
designability1: the minimum designability for one node
designability2: total designability
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include "design.h"
#include "dynamics.h"

int counter, counter2;

void print_help ()
{
	printf ("Usage         : ./checktraj1 mode mode1 mode2 irflag N L inputfile\n");
	printf ("                (./checktraj1 0 0 0 1 11 12 Budding.txt)\n");
	printf ("mode          : 0 normal ; 1 cycled ; 2 triggled cycled\n");
	printf ("mode1         : 0 disallow self-activation ; 1 allow self_activation\n");
	printf ("mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II\n");
	printf ("irflag        : 1 irreducible solution only\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n\n");
	printf ("Output        : activity min_link designability1 designability2\n");
	printf ("activity      : number of 1 in the process\n");
	printf ("min_link      : minimum links of the qualified network\n");
	printf ("designability1: the minimum designability for one node\n");
	printf ("designability2: total designability\n");
}

void init_design (struct designability *design, int argc, char *argv[])
{
	int N, L, i;

	design->mode = atoi (argv[1]);
	design->mode1 = atoi (argv[2]);
	design->mode2 = atoi (argv[3]);
	design->irflag = atoi (argv[4]);
	design->bdflag = 0;
	design->dgflag = 0;
	design->pnflag = 1;
	design->svflag = 1;
	design->N = atoi (argv[5]);
	design->L = atoi (argv[6]);

	design->N += design->mode1;
	N = design->N;
	L = design->L;

	/* the degree of irreducible solution is zero */
	if (design->irflag)
		design->dgflag = 0;

	design->pathway = calloc (sizeof (int), N * L);
	design->NS2 = calloc (sizeof (int), N);
	design->min_link = calloc (sizeof (int), N);
	design->links_dist = calloc (sizeof (int), N * (N + 1));
	design->Links_dist = calloc (sizeof (double), N * N + 1);
	design->degree_dist = calloc (sizeof (int), N * (2 * N + 1));
	design->Degree_dist = calloc (sizeof (double), 2 * N * N + 1);
	design->solution = calloc (sizeof (struct solution), N);

	for (i = 0; i < N; i++)
	{
		init_solution (&(design->solution[i]));
	}
}

void free_design (struct designability *design)
{
	int i;

	free (design->pathway);
	free (design->NS2);
	free (design->min_link);
	free (design->links_dist);
	free (design->Links_dist);
	free (design->degree_dist);
	free (design->Degree_dist);
	free (design->solution);

	for (i = 0; i < design->N; i++)
	{
		free_solution (&(design->solution[i]));
	}
}

void enumerate_solution (struct designability *design, int N, int select[N], int node, int n, struct state state[n], int flux[design->L][n])
{
	int i, j, k, L;

	L = design->L;

	if (node == N)
	{
		counter ++;

		//printf ("Solution %d\n", counter);

		int M;

		M = design->mode1;

		int network[N - M][N - M];

		for (i = M; i < N; i++)
		{
			j = select[i];
			for (k = M; k < N; k++)
			{
				if (design->mode1 == 1 && k == i && design->solution[i].edge[j * N] == 1)
					network[i - M][k - M] = 1;
				else if (design->solution[i].edge[j * N + k] == 0)
					network[i - M][k - M] = 0;
				else if (design->solution[i].edge[j * N + k] == 1)
					network[i - M][k - M] = 1;
				else if (design->solution[i].edge[j * N + k] == -1)
					network[i - M][k - M] = -1;
			}
		}

		reset_state (n, state);
		generate_dynamics (N - M, network, n, state);
		calc_flux (N - M, n, state);

		int p[L - 1];

		for (i = 0; i < L - 1; i++)
		{
			int traj[N - M];

			for (j = M; j < N; j++)
			{
				traj[j] = design->pathway[i * N + j];
			}

			p[i] = pack (N - M, traj);

			flux[i][state[p[i]].flux1] ++;
			/*
			printf ("State");
			for (j = M; j < N; j++)
			{
				printf (" %d", traj[j]);
			}
			printf (" Flux %4d\n", state[p].flux1);
			*/
		}

		if (state[p[L - 2]].flux1 * 2 >= n && state[p[L - 3]].flux1 * 2 >= state[p[L - 2]].flux1 && state[p[L - 4]].flux1 * 2 >= state[p[L - 3]].flux1 && state[p[L - 5]].flux1 * 2 >= state[p[L - 4]].flux1)
			counter2 ++;
	}
	else
	{
		for (i = 0; i < design->solution[node].curr; i++)
		{
			select[node] = i;
			enumerate_solution (design, N, select, node + 1, n, state, flux);
		}
	}
}

void print_pathway (struct designability *design, int n, int flux[design->L][n])
{
	int i, j, p;

	printf ("c %3d %3d %5d %lg\n", design->activity, design->min_links, design->NS2_min, design->NS);

	for (i = 0; i < design->L - 1; i++)
	{
		printf ("State %2d ", i);
		for (j = design->mode1; j < design->N; j++)
		{
			printf ("%d ", design->pathway[i * design->N + j]);
		}
		printf ("Flux ");

		p = 0;
		for (j = 0; j < n; j++)
		{
			p += flux[i][j];
			if (2 * p >= counter)
				break;
		}
		printf ("%4d\n", j);
	}
}

int main (int argc, char *argv[])
{
	struct designability design;

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	init_design (&design, argc, argv);

	FILE *fp;
	char buf[256];
	int i, j, n, N, L, M;

	N = design.N;
	L = design.L;

	M = design.mode1;

	n = 1;
	for (i = M; i < N; i++)
	{
		n *= 2;
	}

	struct state state[n];

	int flux[L - 1][n];

	init_state (N - M, n, state);

	fp = fopen (argv[7], "r");

	while (fgets (buf, 256, fp))
	{
		design.activity = 0;
		for (i = 0; i < L; i++)
		{
			for (j = 0; j < N; j++)
			{
				if (design.mode1 == 1 && j == 0)
				{
					design.pathway[i * N + j] = 1;
					continue;
				}
				fscanf (fp, "%d", &(design.pathway[i * N + j]));
				if (i <= L - 3 || (i == L - 2 && design.mode != 2))
					design.activity += design.pathway[i * N + j];
			}
			fgets (buf, 256, fp);
		}

		if (design.mode == 2)
			design.activity --;

		/* Get the number of solution using enumerating algorithm */
		if (calc_design (&design) == 0)
			continue;

		//print_pathway (&design);

		if (design.svflag && design.pnflag)
		{
			int select[N];

			counter = counter2 = 0;

			for (i = 0; i < L - 1; i++)
			{
				for (j = 0; j < n; j++)
					flux[i][j] = 0;
			}

			enumerate_solution (&design, N, select, design.mode1, n, state, flux);

			print_pathway (&design, n, flux);

			printf ("Irredicible Solution with Good Dynamics Property : %4d / %4d\n", counter2, counter);
		}
	}

	fclose (fp);

	free_design (&design);

	exit (0);
}
